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Abstract We study strategies for increasing the precision in the blurring models by 
maintaining a complexity in the related numerical linear algebra procedures (matrix- 
vector product, linear system solution, computation of eigenvalues etc.) of the same 
order of the celebrated Fast Fourier Transform. The key idea is the choice of a suitable 
functional basis for representing signals and images. Starting from an analysis of the 
spectral decomposition of blurring matrices associated to the antireflective boundary 
conditions introduced in [S. Serra Capizzano, SIAM J. Sci. Comput. 25-3 pp. 1307- 
1325], we extend the model for preserving polynomials of higher degree and fast 
computations also in the nonsymmetric case. 

We apply the proposed model to Tikhonov regularization with smoothing norms 
and the generalized cross validation for choosing the regularization parameter A se- 
lection of numerical experiments shows the effectiveness of the proposed techniques. 

Keywords matrix algebras and fast transforms • Tikhonov regularization • boundary 
conditions 

Mathematics Subject Classification (2000) 65F22 65R32 65T50 



1 Introduction 

We consider the de-convolution problem in the case of signals where the convolution 
kernel is space invariant. In that case the observed signal g : ^ — > M, ^ C M is 
expressible as 
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where / denotes the true signal. The approximation of the integral operator via an 
elementary rectangle formula over an equispaced grid with n nodes leads to a Unear 
system with n equations. 

When imposing proper boundary conditions, the related undetermined linear sys- 
tem becomes square and invertible and fast filter algorithms of Tikhonov type can be 
employed. When talking of fast algorithms, given n the size of the related matrices, 
we mean an algorithm involving a constant number (independent of n) of fast trigono- 
metric transforms (Fourier, sine, cosine, Hartley transforms) so that the overall cost 
is given by 0{nlogn) arithmetic operations. 

For instance, when dealing with periodic boundary conditions, we obtain circu- 
lant matrices which are diagonalizable by using the celebrated fast Fourier transform 
(FFT). Unfortunately such boundary conditions are not always satisfactory from the 
viewpoint of the reconstruction quality. In fact, if the original signal is not periodic 
the presence of ringing effects given by the periodic boundary conditions spoils the 
precision of the reconstruction. 

More accurate models are described by the reflective [11] and antireflective [13] 
boundary conditions, where the continuity of the signal and of its derivative are im- 
posed, respectively. However the fast algorithms are appUcable in this context only 
when symmetric point spread functions (PSFs) are taken into consideration. 

The PSF represents the blur of a single pixel in the original signal. Therefore, 
since it is reasonable to expect that the global light intensity is preserved, the PSF 
is nothing else that a global mask having nonnegative entries and total sum equal 
to 1 (conservation law). Often in several application such a PSF is synometric and 
consequently the symbol associated to its mask is an even function. 

Usually the antireflective boundary conditions lead to better reconstructions since 
linear signals are reconstructed exactly, while the periodic boundary conditions ap- 
proximate badly a linear function by a discontinuous one and the reflective ones by 
a piece-wise linear function: in both the latter case Gibbs phenomena (called ringing 
effects) are observed which are especially pronounced for periodic boundary condi- 
tions. The evidence of such fact is observed in several papers in the literature, e.g. [1, 
3,4,6,12,13]. 

Such good behavior of the antireflective boundary conditions comes directly from 
their definition [13|, since the continuity of the first derivative of the signal was au- 
tomatically imposed. From an algebraic viewpoint, the latter property can be derived 
from the spectral decomposition of the coefficient matrix in the associated linear 
system. Indeed, when considering a symmetric PSF and antireflective boundary con- 
ditions, the hnear system is represented by a matrix whose eigenvalues equal to 1 
(the normalization condition of the PSF coming from the conservation law) are asso- 
ciated to an eigenvector basis spanning all linear functions sampled over a uniform 
grid with n nodes, see [1,3]. In [1], such a remark has been the starting point for 
defining and analyzing the antireflective transform and for designing fast algorithms 
for the spectral filtering of blurred and noisy signals. This algebraic interpretation is 
useful because it can be used for proposing generalizations that preserve the possibil- 
ity of defining fast algorithms, while increasing the expected reconstruction quality 
especially when smooth or piece- wise smooth signals are considered. 
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In this paper, starting from the previous algebraic interpretation, we define higher 
order boundary conditions. This can be obtained by algebraically imposing that the 
spanning of quadratic or cubic polynomials over a proper uniform gridding are eigen- 
vectors related to the normalized eigenvalue 1 . Our proposal improves the antireflec- 
tive model when the true signal is regular enough close to the boundary. Moreover, an 
important property of the proposed approach is that it allows to define fast algorithms 
also in the case of nonsymmetric PSFs (such as the blurring caused by motion). We 
note that reflective and antireflective boundary conditions can resort to fast transforms 
only in the case of symmetric boundary conditions, while in the case of nonsymmet- 
ric PSF we have fast transforms only for periodic boundary conditions which usually 
provide poor restorations for nonperiodic signals. 

In general, if some information on the low frequencies of the signal to be recon- 
structed are available, it is sufficient to impose such sampled components as eigen- 
vectors of the blurring operator related to the eigenvalue 1 (we recall that the global 
spectrum will have 1 as spectral radius). In such a way these component will be main- 
tained exactly by the filtering algorithms since they cut only the spectral components 
related to small eigenvalues (somehow close to zero) which are presumed to be essen- 
tially associated to the noise. In reality, the noise by its random nature of its entries 
will be decomposable essentially in high frequencies while the true signal is sup- 
posed to be approximated in the complementary subspace of low frequencies. There- 
fore, when applying filtering algorithms, if the blurring operator has non-negligible 
eigenvalues associated only to low frequencies (for instance low degree polynomi- 
als), then the reconstruction of the signal will be reasonably good while the noise 
will be efficiently reduced. 

Given this general context, the present note is aimed to define spectral decompo- 
sition of the blurring matrix such that the related transform given by the eigenvectors 
is fast, the conditioning of the transform is moderate (for such an issue in connection 
with the antireflective transform see [5]), and the low frequencies are associated only 
to non-neghgible eigenvalues. 

The organization of the paper is as follows. Section 2 we introduce the deblurring 
problem investigating the spectral decomposition of the coefficient matrix for the 
different kinds of boundary conditions. In Section 3 we define higher order bound- 
ary conditions starting from the spectral decomposition of the antireflective matrix. 
Such transforms are used in Tikhonov-like procedures in Section 4. Section 5 deals 
with a selection of numerical tests on the de-convolution of blurred and noisy signals 
and images. In Section 6 the proposals are extended to a multi-dimensional setting. 
Finally Section 7 is devoted to concluding remarks. 

2 Boundary conditions and associated coefficient matrices 

In this section we introduce the objects of our analysis and we revisit the spectral de- 
composition of blurring matrices in the case of periodic, reflective, and antireflective 
boundary conditions. 

Let f = (. . . ,/o,/i, . . . ,/„,/„+!,. . . )^ be the true signal and {7}"=i the set of in- 
dexes in the field of view. Given a PSF h = {h-m, ...,ho,... ,hm), with 2m + l <n, 
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we can associate to the PSF the symbol 

m 

z{t)= £ V^', i = ^^^. (2.1) 

j=-m 

2. 1 Periodic and Reflective boundary conditions 
Periodic boundary conditions are defined imposing 

f\-j= fn+ 1-j and f„+j = fj , 

for j = l,...,n. The blurring matrix associated to periodic boundary conditions is 
diagonaUzed by the Fourier matrix 

4^ = -^exp(^ n j' ^'^ = l'-'«- 

More precisely, the blurring matrix is 

Ap = (F("))«diag(z(x))FW, (2.2) 

where x; = 2(j — l)7r/n, for ; = l,...,n. We note that the eigenvalues A,- = z(x,) can 
be easy computed by A,- = [F(")(Apei)],/[F(")ei],-, where ei is the first vector of the 
canonical base. 

Reflective boundary conditions are defined imposing 

fi-j = fj and f„+j = fn+i-j, 

for j = If the PSF is symmetric, i.e., h-j = hj, then the blurring matrix 

associated to reflective boundary conditions is diagonaUzed by the cosine transform 
(see [11]) 

where 5,;i = 1 if j = 1 and zero otherwise. More precisely, the blurring matrix is 

A« = (cW)^diag(z(x))cW, (2.3) 

where x, = {i— l)jr/n, for i = I,... ,n. Like for periodic boundary conditions, the 
eigenvalues can be easy computed by A,- = [cW(Aei)],/[cWei],. 
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2.2 Antireflective boundary conditions 

Antireflective boundary conditions are defined imposing (see [13]) 

fi-j = 2/i - fj+i and fn+j = 2f„ - fn-j, 

for; = 

Let Q be the sine transform matrix of order n — 2 with entries 
The antireflective transform of order n can be defined by the matrix (see [1]) 












T = 


P 


Q 


Jp 












where 

_ / «(2n-l) / i-l\ 

for / = 1 , . . . , n and where the permutation matrix / has nontrivial entries Ji^n+i-i = 1. 
i = 1, . . . ,n. We note that ||p||2 = 1; moreover/ is often cafled flip matrix. 

If the PSF is symmetric and 2w + 1 < n — 2, the spectral decomposition of the 
coefficient matrix in the case of antireflective boundary conditions is 

AA = Tdmg{z{y))T-\ (2.5) 

with y defined &syi = {i — i)n/{n — I) fox i= I,..., n — I andjn = 0. The eigenval- 
ues of A can be computed in O(nlogn) real operations resorting to the discrete sine 
transform (see [2]). 

Concerning the inverse antireflective transform in [1] we have given its ex- 
pression and the resulting form is analogous to that of the direct transform T. As a 
matter of fact, given an algorithm for the direct transform, a procedure for computing 
the inverse transform needs only to have a fast way for multiply J"' by a vector. 

Remark 2.1 Observe that ,5^i = span{p,/p} is the subspace spanned by equi-spaced 
samplings of linear functions. In that case its linear complement is given by S^f- = 
span{sin(7x)}"~j, with Xi = {i — l)n/{n — I), i = I, . . . ,n. Unfortunately such a Un- 
ear complement is not orthogonal and consequently the related transform cannot be 
unitary, as long as we maintain such a trigonometric basis useful for the fast com- 
putations. Up to standard normalization factors this choice leads to the antireflective 
transform (2.4). 

Remark 2.2 Implementing filtering methods, Uke Tikhonov, J^/ is about fuUy pre- 
served since the associated eigenvalues are z(0) = 1. 
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3 Higher order boundary conditions 

Starting from Remarks 2.1 and 2.2, we define higher order boundary conditions 
which represent the main contribution of this work. The approach in Section 2 de- 
fines accurate boundary conditions imposing a prescribed regularity to the true signal 
/. The study of the spectral decomposition of the associated coefficient matrices is 
a subsequent step for defining fast and stable filtering methods. In this section, we 
define higher order boundary conditions starting from the eigenspace, i.e., the signal 
components, that we wish to preserve. 

We start by imposing to preserve S^i and by suggesting other choices for S^f . 
By the way, the request of giving fast algorithms suggests the use of a cosine or 
exponential basis in place of that of sine functions, both for the direct and inverse 
transforms. To preserve polynomials of low degree and at the same time to resort 
to fast trigonometric transforms, we need a transform with a structure analogous to 
(2.4). Therefore we need the cosine transform and the Fourier matrix of order n — 2. 
We define F = F^"-^) and C = &'-^\ explicitly 



C, = ./^cos^('--l)(^^"-l)^ 



In -A 



Fij = —i== exp 



and 

^^___'-i27r(i-l)0--l) 

for /,./'= 1, ... ,n — 2. 

We note that the first column of and of F^ are a samphng of the constant 
function. Hence the span of the columns of or F^ has a nontrivial intersection 
with 5^1. Accordingly, we choose the two vectors for completing these trigonomet- 
ric basis as a uniform sampling of a quadratic function instead of a linear function. 
More precisely, instead of S^i we consider 5^q = span{q,7q}, where q is a uniform 
sampling of a quadratic function in an interval that will be fixed later 

The interval and the samphng grid for the basis functions of our transform are 
fixed according to the following remark. 

Remark 3.1 Up to normalization, the y'th column of Q is sin(yx), where x, = in/{n — 
1), for i = 1, ... ,n — 2. Extending the sampling grid such that the jth frequency is 
extended by continuity, we add the grid points xq = and i = K. Since sin(jxo) = 
sin(7x„_ 1 ) = for all j, we obtain exactly the two zero vectors in the first and the last 
row of T , i.e., the (/' + l)th column of T is the jth column of Q extended in xq and 
Xn, for y = 1 , . . . , n — 2. The first and the last column of T are the sampling of hnear 
functions at the same equispaced points x,- e [0, tt] , i = 0, . . . , n — 1 . 



3.1 The case of symmetric PSF 

Firstly, we consider a symmetric PSF. In such case we can use the cosine basis. Up to 
normaUzation, the (7-|-l)thcolunnnof is cos(7x), where x,- = (2j— l);r/(2n — 4), 
for J = 1,. . . ,n — 2. Extending the grid by continuity, we addxo = — ;r/(2n — 4) and 
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x„_ 1 (2n — 3 ) tt/ (2n — 4) . With this extended grid we can define the basis functions 
as a n points uniform sampUng of the interval 



[a, b] = 



7t (2n-3)n 



In- A' 2n-4 



where the grid points are 

Xi = {2i - l)n/{2n - 4), i = 0,...,n-l. 



(3.1) 



(3.2) 



We fix q = q/||q||2, where [q],+i = (b — Xi)^, i = 0,...,n—l. The fast transform 



associated to and can be defined as follows: 







T 

C 




Tc = 


q 




Jq 






T 





(3.3) 



with [Ca]j = y ^;^cos((j - l)fl) and [c^,]; - y -;p2-^"^vu- - V- v waij 
since ^ = TT — a, for j = 1 — 2. 

It remains to define the eigenvalues associated to Tc. Since we want to preserve 
S^q, similarly to what was done for p in the case of the antireflective boundary con- 
ditions, we associate to q and Jq the eigenvalue z(0) = 1. Concerning the other fre- 
quencies, since they are defined by the cosine transform, we consider the eigenvalues 
of the reflective matrix in (2.3), but of order n — 2. 

In conclusion, for the case of a symmetric PSF, we define a new blurring matrix 
using the following spectral decomposition 



V^cos((y-i)fc) = (-i)^-i[c«] 



Ac 



Tcdmg{z{x))T-\ 



(3.4) 



where x,- = {i — 2)n/{n — 2), for j = 2, . . . ,n — 1, and xi = x„ = 0. 

We note that z{xi) ~ z(x„) = 1, while the eigenvalues for = 2, . . . ,n — 1, 
are the same of Ar of order n — 2 and hence they can be computed in 0(«Iogn) by 
a discrete cosine transform. The product of Tc by a vector can be computed mainly 
resorting to the inverse discrete cosine transform. The inverse of Tc will be studied 
in Subsection 3.3, where we will show that the product of by a vector can be 
computed mainly resorting to a discrete cosine transform. Therefore the spectral de- 
composition (3.4) can be used to define fast filtering methods in the case of symmetric 
PSFs. Moreover, we expect an improved restoration with respect to the antireflective 
model since Ac preserves uniform sampUngs of quadratic functions while Ar pre- 
serves only uniform samplings of linear functions. 



3.2 The case of nonsymmetric PSF 

In the case of nonsyrmnetric PSF we can use the exponential basis. Up to normaliza- 
tion, the (7-|-l)th column of is exp(i7x), where i = \/^andx,- = {i—l)2n/{n — 
2), for i = 1, ... ,n — 2. Extending the grid by continuity, we add xq = —2n/{n — 2) 
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and x„_ i=2n. With this extended grid, we can define the basis functions as a n points 
uniform sampUng of the interval 

[a,b] = [-2K/{n-2),27t], (3.5) 

where the grid points are 

Xi = {i-l)2n/{n-2), j = 0,...,n- 1. (3.6) 

We note that the interval and the grid points in the nonsymmetric case are dif- 
ferent with respect to the symmetric case (compare (3.5) with (3.1) and (3.6) with 
(3.2)). Therefore, defining q = q/||q||2, where [q],+i = (b — Xi)^, i = 0,. . . ,n— I, it is 
different from which obtained in the symmetric case in the previous subsection. The 
fast transform associated to and can be defined as follows 



(3.7) 



where [Ca]j = exp{i{j— l)a)/\/n — 2and [c^jy = exp(i(7 — l)fo)/\/n^-2= l/\/n — 2, 
for^' = l,...,n — 2. 

It remains to define the eigenvalues associated to 7>. Similarly to what done for 
Tc, we associate to q and Jq the eigenvalue z(0) = 1, while for the other frequencies 
we consider the eigenvalues of the circulant matrix in (2.2), but of order n — 2. 

Consequently, in the case of a generic PSF, we define a new blurring matrix using 
the following spectral decomposition 

= 7>diag(z(x))r/i, (3.8) 

where a:,- = {i — 2)2n/n, for i = 2, . . . ,n — 1, andxi =x„=0. 

We note that =z{x„) = 1, while the eigenvalues z(x,), for / = 2, . . . ,n — 1, are 
the same of Ap of order n — 2 and hence they can be computed in O(nlogn) by a fast 
Fourier transform. The product of Tp by a vector can be computed essentially resort- 
ing to the inverse fast Fourier transform. The inverse of 7> will be studied together 
with the inverse of Tc in the next Subsection, where we will show that the product of 
Tf^ by a vector can be computed by using the fast Fourier transform. Therefore the 
spectral decomposition (3.8) can be used to define fast filtering methods also in the 
case of nonsymmetric PSFs. 



3.3 The inverse transform 

In this subsection we show that the inverse of Tc and the inverse of Tp are fast trans- 
forms. This means that the associated matrix vector product can be performed mainly 
via a suitable trigonometric transform. 











Tf = 


q 










H 





Fast transforms for high order boundary conditions 



9 



Theorem 3.1 Let 







T 

C 




Tx = 


q 















(3.9) 



be a given nxn matrix, where J is the flip matrix and X is a discrete trigonometric 
transform such that JXJ = X. Then T^^y can be computed in 0{nlog{n)) for all 
yeC. 



Proof We note that 



where 



Tx = Tx + [ei\en] 



Oci 
cl 0, 



(3.10) 







0^ 




Tx = 


q 




Jq 






0^ 





is easy to invert. Hence T^^ can be computed by the Sherman-Morrison- Woodbury 
formula. 

We compute Since qn = the first and the last row can be decoupled and 
we look for of the form 





a 


0^ 


" 


1 _ 


V 


X 


J\ 







0^ 


a 



Tx 



Fixing q = [^1 , q^, 0]^, by direct computation a = \/qi and v = —Xq/q\. Therefore, 
V can be computed in O(nlogn) by a trigonometric transform. For the implementation 
it can be explicitly computed and inserted into the code. 

Given A e C"^" and i7,y e C"^*^, the Sherman-Morrison- Woodbury formula is 

[9]: 

{A + UV")~^ =A-^ -A-^U{I + V"A-^U)~^V"A-^. (3.11) 

It can be very useful for computing the inverse of A -|- UV^ when kin, taking into 
account the possible instability. Applying the formula (3.1 1) to (3.10) we obtain 



a 
v J\ 
a 



1 + 



cl 



Tx [ei|e„] 



-1 



0c„^0 
cl 



c^v l+clJ\ 



n -1 



We note that c^X and clX can be computed in 0(nlog(n)) and moreover they can 
be explicitly computed and inserted into the implementation like done for the vector 
v. In this way the matrix vector product for ' requires a fast discrete trigonometric 
transform of 0(nlog(n)) plus few lower order operations between vectors. □ 

From Theorem 3.1, it follows that the product of T^^ and by a vector can 
be computed in 0(nlog(n)) and hence they are fast transforms. 
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4 Tikhonov regularization with fast transforms 

We consider the Tikhonov regularization, where the regularized solution is computed 
as the solution of the following minimization problem 



(4.1) 



where, jU is the properly chosen regularization parameter, g is the observed signal, A 
is the coefficient matrix and L is a matrix such that Null (A ) f) Null (L) = (see [7]). 
The matrix L is usually the identity matrix or an approximation of partial derivatives. 

It is convenient to define L using the same boundary conditions of A in order 
to obtain fast algorithms. For instance, L equal to the Laplacian with antireflective 
boundary conditions is 

^0 ... 
-1 2 -1 







-1 2 -1 




(4.2) 



: because Null(L) = 
^ for s{x) = (2 — 2cos(x)) and y defined accord- 

-.(0)=0). 



We note thatdim(Null(L)) = 2. However, Null(A)nNull(L) 
Yi. We have La = Tx dmg{s{y) ) 
ing to (2.5) (note that s{yi ) .?(>'„) 

Using the approach in Section 3, for high order boundary conditions the Laplacian 
matrix can be defined similarly by 



Lc = Tcdiiig{s{y))Tc\ 
Lf =7>diag(i(x))7>-i, 



where the grid points y and x are defined according to (3.4) and (3.8) respectively. 



4.1 Tikhonov regularization and reblurring 

The minimization problem (4.1) is equivalent to the normal equations 

(A^A + /iL^L)f = A^g. (4.3) 

Regarding the antireflective algebra, in [6] it was observed that the transposition op- 
eration desttoys the algebra structure and leads to worse restorations with respect to 
reflective boundary conditions. To overcome this problem, in [4] the authors proposed 
the reblurring which replaces the transposition with the correlation operation. More- 
over, it was shown that the latter is equivalent to compute the solution of a discrete 
problem obtained by a proper discretization of a continuous regularized problem. 
Practically, we replace A^ with A' obtained imposing the same boundary conditions 
to the coefficient matrix arising from the PSF rotated by 180 degrees. The matrices 
defined in (2.4), (3.4), and (3.8) can be denoted by Ax{z), X e {A,C,F} since they 
are univocaUy defined by the function z when the transform Tx is fixed. With this 
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notation (z) ^Ax(z) since the rotation of the PSF exchange hj with in (2.1), 
which corresponds to take z. Therefore, in the case of periodic boundary conditions 
A' = A^, but this is not true in general for the other boundary conditions. If the PSF 
is symmetric then A' =A for every boundary conditions. 

In the following we use the reblurring approach and hence we replace (4.3) with 

(A'A + ML'L)f,eg=A'g. (4.4) 

If we use the same boundary conditions for A and L, the (4.4) can be written as 
+ Ml*P)freg =Ax{z)g. In [5] it is proved that for antireflective boundary con- 
ditions (4.4) defines a regularization method when L = I and the PSF is symmetric. 

If the spectral decomposition of A is A = TxDT^^, where D = diag(d), and L = 
Txdiag{s)T^^ , then the spectral filter solution in (4.4) is given by 

f.eg = Tx0D-%-^ g, = diag,..,...,„ (^^p^) (4.5) 



4.2 GCV for the estimation of the regularization parameter 

A largely used method for estimating the regularization parameter is the general- 
ized cross validation (GCV) [8]. For the method in (4.5), GCV determines the regu- 
larizing parameter ^ that minimizes the GCV functional 

^ ' trace(/-A73f0£»-ir-i)2 

where freg is defined in (4.5). For Ap (periodic boundary conditions) and Ar (re- 
flective boundary conditions with a symmetric PSF), the equation (4.4) is exactly 
equation (4.3). In this case, in [10] it is proven that 

G(;i) = ^kMf, (4.7) 

where a,- = \si\^/{\di\^ + ju|i,f ) and g = T£^g, for equal to and CW, re- 
spectively. 

Here we have 

||g-Afreg||2= ||7i(/-^)r^-'g||2 

« ||(/-0)r^-ig||2, (4.8) 

because Tx is not unitary but it is "close" to a unitary matrix since it is a rank four 
correction of a unitary matrix. For the estimation of the SVD of T (antireflective 
boundary conditions case) see [5]. Therefore, we compute the regularization param- 
eter ^ by minimizing the same functional as in (4.7). More precisely 
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0.2 0.4 0.6 O.B 1 



Fig. 5.1 tme signal, — observed signal with Gaussian blur and 0.1% of noise. The vertical lines 

denote the field of view. 

5 Numerical experiments 

We present some signal deblurring problems. The restorations are obtained by em- 
ploying Tikhonov regularization using (4.5) with smoothing operator L equal to the 
Laplacian. The code is implemented in Matalab 7.0. 

In the first example the observed signal is affected by a Gaussian blur and 0.1% 
of Gaussian noise. True and observed signals are shown in Figure 5.1. We consider 
a low level of noise because in such case the restoration error is mainly due to the 
error of the boundary conditions model. Since the PSF is symmetric, we compare our 
blurring matrix Ac with reflective and antireflective boundary conditions. 

Let f be the true signal, the relative restoration errors (RRE) ||f — freglb/ljflU is 
plotted in Figure 5.2. In such figure it is evident that Ac provides restorations with 
a lower RRE with respect to antireflective boundary conditions, which are already 
known to be more precise than reflective boundary conditions. Moreover, the RRE 
curve varying the regularization parameter /i is flatter with respect to the other bound- 
ary conditions. This allows a better estimation of the regularization parameter using 
the GCV. The value /iocv that gives the minimum of the GCV functional in (4.7) 
is reported in Figure 5.2 with a It is evident that in the case of Ac the RRE ob- 
tained with Hgcm is closer to the minimum with respect to antireflective boundary 
conditions. The minimum RRE is 0.135 for Ac while it is 0.177 for the antireflective 
boundary conditions. Moreover, for Ac we obtain /Igcv = 5.6 x 10^^ which gives a 
RRE equal to 0.135, while for antireflective boundary conditions jUgcv = 1.6x lO^'* 
gives a RRE equal to 0.502. 

The quality of the restoration is validated also from the visual evidence of the 
restored signals. In Figure 5.3 we show the restored signal corresponding to /iopt, 
which is the value of the regularization parameter pL corresponding to the minimum 
RRE, and to /iocv ■ We note that Ac gives a better restoration especially for preserving 
jumps in the signal. On the other hand this implies a slightly lose of the smoothness 
of the restored signal. Eventually, using /Xqcv our proposal with Ac gives a good 
enough restoration while this is not true for the antireflective boundary conditions. 

The second example is a moving PSF with a 1% of Gaussian noise. True and ob- 
served signals are shown in Figure 5.4. Since the PSF is nonsymmetric we consider 
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10^ 10" 10"^ 10° 



Fig. 5.2 RRE: — Ac, antireflective, - ■ - reflective (* denotes values coiTesponding to ^iocv)- 



1.2 




'0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Fig. 5.3 Restored signals: — Ac, ■ ■ ■ antireflective, true signal. 




Fig. 5.4 — true signal, observed signal with moving blur and 1% of noise. The vertical lines denote 

the field of view. 

Af instead of Ac- Moreover, since antireflective and reflective boundary conditions 
lead to matrices that can not be diagonalized by fast transforms, we can compare 
Ap only with periodic boundary conditions. From Figure 5.5 and 5.6 we note that 
the same considerations done in the previous example hold unchanged. Indeed the 
minimum RRE is 0.091 for A/r while it equals 0.198 for the periodic boundary con- 
ditions. Moreover, for A/r we obtain Hocw = 1-7 x 10^^ which gives a RRE equal to 
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Fig. S.S RRE: — Af, periodic (* denotes values con'esponding to Hacy)- 




a 0.2 




0.6 0. 




0.2 0. 



Mopt 

Fig. S.6 Restored signals: — Af, ■ ■ ■ periodic, true signal. 



Mgcv 



0.096, while for periodic boundary conditions jUgcv = 2.7 x 10 gives a RRE equal 
to 0.704. 



6 The multidmensional case 



A standard way for defining the multidimensional transform is by tensor product. 
Thus 



Ad) 



rid) 



)73f. 



d times, where n= (ni, . . . and Tx.m is the transform Tx of order m. For a 2D 
array of size n x m, this is easily implemented doing m ID transforms of size n for 
each column and then n ID transforms of size m for each row. 

The computation of the eigenvalues is more involved. The strategy is the same 
described in [2] for computing the eigenvalues of antireflective matrices. The algo- 
rithm in Section 3.2. 1 in [2] can be applied to our proposal, by replacing the discrete 
sine transform by the cosine or the Fourier transform. More in details, in the 2D case: 

1. Compute two ID PSF summing the rows and the columns of the 2D PSF. 

2. Apply two ID transforms Tx for computing the eigenvalues that correspond to the 
frequencies indexed as the edges of the image (the vertical edges are associated 
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(a) (b) 

Fig. 6.1 (a) True image, (b) Observed image with out of focus blurring and 0.1% of Gaussian noise. 



10" 




10^ 10"^ 



Fig. 6.2 RRE: — Ac, antireflective, - ■ - reflective (* denotes values corresponding to /iocv)- 

to the PSF obtained summing the columns and the horizontal edges to the other 
PSF). 

3. Apply a 2D cosine or Fourier transform for computing the eigenvalues indexed 
as the inner part of the image. 



6.1 Image deblurring 

We consider the deblurring problem with an out of focus blur and 0.1% of Gaussian 
noise. The true and the observed images are shown in Figure 6.1. The restored images 
are obtained by using the smoothing operator L = I. 

We note that Ac gives a better restoration with respect to antireflective and reflec- 
tive boundary conditions (see Figure 6.2). In Table 6.1 the RRE is shown for jUopt and 
Mgcv. while in Figures 6.3 and 6.4 we have the restored images for the considered 
boundary conditions and the two choices of jj.. 

Even if there is not a large reduction of the RRE, the images restored with Ac 
show lesser ringing effects with respect to the antireflective boundary conditions at 
least in the south-west corner of the image. 
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Table 6.1 RRE for the restoration of the observed image in Figure 6.1. 







Mgcv 


reflective 


0.0647 


0.0723 


antireflective 


0.0570 


0.0602 


Ac 


0.0564 


0.0579 




(a) reflective (b) antireflective (c) Ac 

Fig. 6.3 Restored images for n^^. 




(a) reflective (b) antireflective (c) Ac 

Fig. 6.4 Restored images for ^tccv- 

For a general image the use of Ac instead of antireflective boundary conditions 
leads to negligible improvement if the image is not smooth enough at the boundary 

or if the noise level is so high to dominate the approximation error in the restoration. 

For concluding, we consider a nonsymmetric PSF. The observed image in Figure 
6.5 (a) is affected from an out of focus combined with a moving blur. Since the PSF 
is nonsymmetric, we compare with periodic boundary conditions like in Section 
5. In Figure 6.5 (b) the RRE for Ap is significantly lower than the RRE of periodic 
boundary conditions. Indeed, in Figure 6.6 it is possible to note that in the case of 
periodic boundary conditions, the ringing effects at the edges (in the direction of the 
motion) damage completely the restoration also for /iopt- Moreover, the GCV gives a 
good estimation of the regularization parameter only in the case of as it is evident 
in the plot of Figure 6.5 (b). 
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10° 




7 Conclusions 



In Section 3 we have given a framework to construct precise models for deconvolution 
problems using fast trigonometric transforms. The same idea could be applied to 
different problems having a shift invariant kernel. Indeed, if we have information on 
the signal to restore, the set S^i can be replaced by other functional spaces that we 
want to preserve. Moreover, higher order boundary conditions can be constructed, 
even if the numerical results show that for image deblurring problems this approach 
does not give substantial improvements. 

The introduced fast transforms was applied in connection with Tikhonov regu- 
larization and the reblurring approach. However, they could be useful also for more 
sophisticated regularization methods like Total Variation for instance. 

The analysis of the Tikhonov regularization in Section 4 is useful also for the 
antireflective boundary conditions. Indeed, it was not previously considered in the 
literature the case ofL^I and the choice of the regularization parameter /i using the 
GCV. 

Since the proposed transforms are not orthogonal, they were applied in connec- 
tion with the reblurring approach, but the theoretical analysis of the regularizing prop- 



18 



Marco Donatelli 



erties of such approach exists only in the case of antireflective boundary conditions 
and symmetric kernel (see [5]). Therefore, a more detailed analysis, especially in 
the multidimensional case with a nonsymmetric kernel, should be considered in the 
future. 

Acknowledgements I would thank Serra Capizzano for useful discussions. 
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